function [parms_fitted, gesol, FitVal, FLAG_SUCCESSFUL_FSOLVE] = ...
    FitModel20220926_ver3(parms, kappa_guess, b_guess, eta_guess, z1_guess, a_guess, b_sigma_x_REC, ...
        U_mean_target, U_std_target, wagebillratio_target, vol_gy_target, ave_mu3_vol_target, NumTries, fsolve_opts)
% target moments:
% E[U], vol[u], vol(N*w)/vol(Y), vol(dlogY), vol(ylds), vol(mu3_1y)
% parameters:
% kappa, b, eta, z(1), chi, sigma_x_REC
% sigma_x(N,z) = exp(a + sigma_x_REC*1{z==1})

    if z1_guess>=0; error('z1_guess should be <0'); end
    if eta_guess<=0 || eta_guess>=1; error('check eta_guess'); end
    if kappa_guess<=0; error('check kappa0_guess'); end
    if b_guess<=0; error('check b_guess'); end

    if parms.Nz~=2; error('Code is for 2 states.'); end
    if abs(parms.grid_z(2)-0)>1e-10; error('Code assumes z(2)=0.'); end
    
%     fixedpt_tol = 1e-8;
    % fixedpt_tol = 1e-7;
    fixedpt_tol = 1e-6; % speed things up
    
    N_Nxz = repmat(parms.grid_N(:),[1 parms.Nx parms.Nz]);
    N_Nz = repmat(parms.grid_N(:),[1 parms.Nz]);
    U_Nz = 1 - N_Nz;
    U_Nxz = 1 - N_Nxz;
    x_Nxz = permute(repmat(parms.grid_x(:),[1 parms.NN parms.Nz]),[2 1 3]);
    
    phi_of_x = @(x) 1./(1 + exp(-x));
    phi_Nxz = phi_of_x(x_Nxz);
    log_phi_Nxz = log(phi_Nxz);
    input_mu3_Nxz = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    
    FLAG_SUCCESSFUL_FSOLVE = false;
    
    for iter = 1:NumTries
        fprintf('iter %g with kappa_guess=%g, b_guess=%g, eta_guess=%g, z1_guess=%g, a_guess=%g.\n', ...
            iter, kappa_guess, b_guess, eta_guess, z1_guess, a_guess);
        [XFit, FitVal, EXITFLAG] = fsolve(@FitFun,[log(kappa_guess); ...
            log(b_guess); log(1/eta_guess-1); log(-z1_guess); ...
            a_guess],fsolve_opts);
        if EXITFLAG>0
            fprintf('Found roots on iteration %g.\n', iter);
            FLAG_SUCCESSFUL_FSOLVE = true;
            break;
        else
            fprintf('iter %g, fitted vals (not solution):\n kappa=%g, b=%g, eta=%g, z1=%g, a=%g.\n', ...
                iter, exp(XFit(1)), exp(XFit(2)), 1/(1 + exp(XFit(3))), ...
                -exp(XFit(4)), XFit(5));
            kappa_guess = 0.1 + 5*rand;
            b_guess = 0.3 + 0.6*rand;
            eta_guess = 0.1*rand;
            z1_guess = fit_std_dz_2state(parms,vol_gy_target/sqrt(2.5+1*rand));
            a_guess = log(0.01 + 0.2*rand);
        end
    end
    
    % do a last compute
    parms_fitted = parms;
    parms_fitted.kappa(:) = exp(XFit(1));
    parms_fitted.b(:) = exp(XFit(2));
    parms_fitted.eta(:) = 1/(1 + exp(XFit(3)));
    parms_fitted.grid_z(1) = -exp(XFit(4));
    parms_fitted.sigma_x_a = XFit(5);
    parms_fitted.sigma_x_REC = b_sigma_x_REC;
    parms_fitted.sigma_x(:,1) = exp(XFit(5) + b_sigma_x_REC);
    parms_fitted.sigma_x(:,2) = exp(XFit(5));

    fprintf('kappa=%g, b=%g, eta=%g, z1=%g, a=%g.\n', ...
        parms_fitted.kappa(1), parms_fitted.b(1), parms_fitted.eta(1), ...
        parms_fitted.grid_z(1), XFit(5));

    gesol = SolveFixedPoint_Model_20220926(parms_fitted,fixedpt_tol);
    parms_fitted.FitVal = FitVal;

    if ~FLAG_SUCCESSFUL_FSOLVE
        warning('Model fit: not all moments are satisfied.');
    end

    function Z = FitFun(X)
        
        Z = NaN*zeros(5,1);
        
        parms_temp = parms;
        parms_temp.kappa(:) = exp(X(1));
        parms_temp.b(:) = exp(X(2));
        parms_temp.eta(:) = 1/(1 + exp(X(3)));
        parms_temp.grid_z(1) = -exp(X(4));
        parms_temp.sigma_x(:,1) = exp(X(5) + b_sigma_x_REC);
        parms_temp.sigma_x(:,2) = exp(X(5));
        
        z_Nxz = permute(repmat(parms_temp.grid_z(:),[1 parms.NN parms.Nx]),[2 3 1]); % note: this must be updated when z changes!!
        Y_Nxz = exp(z_Nxz).*N_Nxz;

        gesol_temp = SolveFixedPoint_Model_20220926(parms_temp,fixedpt_tol);
        if gesol_temp.VFIerr>fixedpt_tol
            warning('FitModel20220926_ver1: VFI error = %g is greater than tolerance of %g.', ...
                gesol_temp.VFIerr, fixedpt_tol)
            bComputeMoments = false;
        else
            bComputeMoments = true;
        end
        
        %% Newer version
        
        if bComputeMoments

        [Phi_Nxz,ss_dist_Nxz] = MakeSparseTransMatrix_Nxz(parms_temp,gesol_temp);

        % NEW CODE: quarterly average
        [mean_U, std_U] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,U_Nxz);
        [mean_Y, std_Y] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,Y_Nxz);
        [mean_Nw, std_Nw] = Calc_quarterly_ave_NHz(Phi_Nxz,ss_dist_Nxz,N_Nxz.*gesol_temp.wages);
        ratio_stdNw_stdY = std_Nw/std_Y;
        
        % quarterly output vol (exact), from library\exact moments
        vol_gy = moment_quarterly_output_growth_vol_Nxz(parms_temp,gesol_temp,Phi_Nxz,ss_dist_Nxz);
        
%         gesol_temp = PriceRealBonds_Model20220926(parms_temp,gesol_temp);
        
%         % real yield curve
%         mats = 12:12:60;
%         real_yld_mean = zeros(size(mats));
%         real_yld_std = zeros(size(mats));
%         for n = 1:numel(mats)
%             idx = find(gesol_temp.real_bonds.maturities==mats(n));
%             mat_yr = mats(n)/12;
%             [real_yld_mean(n),real_yld_std(n)] = Calc_monthly_vol(ss_dist_Nxz,-log(gesol_temp.real_bonds.prices(:,:,:,idx))/mat_yr);
%         end
%         real_yld_std = 100*real_yld_std;
        
        std_mu3_cons_1y = CalcTimeSeriesDiffVol(input_mu3_Nxz,ss_dist_Nxz,Phi_Nxz,12);
        
        Z(1) = (mean_U - U_mean_target)/U_mean_target;
        Z(2) = (std_U - U_std_target)/U_std_target;
        Z(3) = (ratio_stdNw_stdY - wagebillratio_target)/wagebillratio_target;
        Z(4) = (vol_gy - vol_gy_target)/vol_gy_target;
        Z(5) = (std_mu3_cons_1y - ave_mu3_vol_target)/ave_mu3_vol_target;
        
        end
        
    end

end